#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### Identity Politics and Trade Preferences #####
#
# Review of International Political Economy
#
# Tyler Girard, Andrea Lawlor, and Erin Hannah
#
#
#
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#

#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### Libraries #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#

library(here)
library(tidyverse)
library(table1)
library(flextable)
library(modEvA)
library(optiscale)
library(circular)
library(fastDummies)
library(lmtest)
library(modelsummary)
library(marginaleffects)
library(ggsignif)


#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### Set Seed #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#

set.seed(515)

#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### Helper Functions #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#

source(here("Functions", "mdpref.R"))

tidy_custom.polr <- function(x, ...) {
  s <- lmtest::coeftest(x)
  out <- data.frame(
    term = row.names(s),
    p.value = s[, "Pr(>|t|)"])
  out
}

#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### Data #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#

df <- rio::import(here("Data", "Trade Policy_June 9, 2023_05.51.csv"))


head(df)
names(df)

df <- df[-c(1,2),]


#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### Test IDs and Multiple Submissions #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#

testPID <- c("c1137edb-8aee-ed11-8e8b-0022486e087a",
             "c3137edb-8aee-ed11-8e8b-0022486e087a",
             "ba137edb-8aee-ed11-8e8b-0022486e087a",
             "bc137edb-8aee-ed11-8e8b-0022486e087a")

dups <- df %>%
  janitor::get_dupes(PID) %>%
  group_by(PID) %>%
  slice(1) %>%
  pull(PID)


#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### Straightliners #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#

#### ~ News ~ ####
straights_news <- df %>%
  filter(!PID %in% testPID) %>%
  filter(!PID %in% dups) %>%
  select(PID, contains("news_")) %>%
  select(-c(news_DO))

head(straights_news)

levels=unique(do.call(c,straights_news[,-1])) #all unique values in df
out <- sapply(levels,function(x)rowSums(straights_news[,-1]==x)) #count occurrences of x in each row
colnames(out) <- levels
out <- as.data.frame(cbind(straights_news, out))

straights_news_scale <- straights_news %>%
  pivot_longer(cols = -c(PID)) %>%
  mutate(value = na_if(value, "-99")) %>%
  mutate(value = case_when(value == "Never" ~ 1,
                           value == "A few times a month" ~ 2,
                           value == "A few times a week" ~ 3,
                           value == "Daily" ~ 4)) %>%
  pivot_wider(names_from = name, values_from = value)

straights_news_scale$newconscale <- apply(straights_news_scale[,-1], 1, function(x) mean(x, na.rm = TRUE))
straights_news_scale$newconscale01 <- modEvA::range01(straights_news_scale$newconscale, na.rm = TRUE)

straights_news_scale <- straights_news_scale %>%
  select(PID, newconscale01)

straights_news_pid <- out %>%
  filter(`-99` == 7 |
           `A few times a month` == 7 |
           `A few times a week` == 7 |
           `Never` == 7 |
           `Daily` == 7) %>%
  pull(PID)

#### ~ Racial Resentment ~ ####
straights_rr <- df %>%
  filter(!PID %in% testPID) %>%
  filter(!PID %in% dups) %>%
  select(PID, contains("racialresent")) %>%
  select(-c(racialresent_DO))

head(straights_rr)

levels=unique(do.call(c,straights_rr[,-1])) #all unique values in df
out <- sapply(levels,function(x)rowSums(straights_rr[,-1]==x)) #count occurrences of x in each row
colnames(out) <- levels
out <- as.data.frame(cbind(straights_rr, out))

straights_rr_pid <- out %>%
  filter(`-99` == 6 |
           `Strongly disagree` == 6 |
           `Somewhat disagree` == 6 |
           `Neither agree nor disagree` == 6 |
           `Somewhat agree` == 6 |
           `Strongly agree` == 6) %>%
  pull(PID)

#### ~ Vulnerability ~ ####
straights_vuln <- df %>%
  filter(!PID %in% testPID) %>%
  filter(!PID %in% dups) %>%
  select(PID, contains("empfin")) %>%
  select(-c(empfin_DO))

head(straights_vuln)

levels=unique(do.call(c,straights_vuln[,-1])) #all unique values in df
out <- sapply(levels,function(x)rowSums(straights_vuln[,-1]==x)) #count occurrences of x in each row
colnames(out) <- levels
out <- as.data.frame(cbind(straights_vuln, out))

straights_vuln_pid <- out %>%
  filter(`-99` == 4 |
           `Strongly disagree` == 4 |
           `Somewhat disagree` == 4 |
           `Neither agree nor disagree` == 4 |
           `Somewhat agree` == 4 |
           `Strongly agree` == 4) %>%
  pull(PID)


#### ~ Exp 1 Matrix 1 ~ ####
straights_expmat1 <- df %>%
  filter(!PID %in% testPID) %>%
  filter(!PID %in% dups) %>%
  select(PID, contains("exp1_mat1")) %>%
  select(-c(exp1_mat1_DO))

head(straights_expmat1)

levels=unique(do.call(c,straights_expmat1[,-1])) #all unique values in df
out <- sapply(levels,function(x)rowSums(straights_expmat1[,-1]==x)) #count occurrences of x in each row
colnames(out) <- levels
out <- as.data.frame(cbind(straights_expmat1, out))

straights_expmat1_pid <- out %>%
  filter(
    `Strongly disagree` == 3 |
      `Somewhat disagree` == 3 |
      `Neither agree nor disagree` == 3 |
      `Somewhat agree` == 3 |
      `Strongly agree` == 3) %>%
  pull(PID)

#### ~ Exp 1 Matrix 2 ~ ####
straights_expmat2 <- df %>%
  filter(!PID %in% testPID) %>%
  filter(!PID %in% dups) %>%
  select(PID, contains("exp1_mat2")) %>%
  select(-c(exp1_mat2_DO))

head(straights_expmat2)

levels=unique(do.call(c,straights_expmat2[,-1])) #all unique values in df
out <- sapply(levels,function(x)rowSums(straights_expmat2[,-1]==x)) #count occurrences of x in each row
colnames(out) <- levels
out <- as.data.frame(cbind(straights_expmat2, out))

straights_expmat2_pid <- out %>%
  filter(`-99` == 3 |
           `Strongly disagree` == 3 |
           `Somewhat disagree` == 3 |
           `Neither agree nor disagree` == 3 |
           `Somewhat agree` == 3 |
           `Strongly agree` == 3) %>%
  pull(PID)

#### ~ All Straights ~ ####

allstraights <- c(straights_news_pid,
                  straights_rr_pid,
                  straights_vuln_pid,
                  straights_expmat1_pid,
                  straights_expmat2_pid)

allstraightsPID <- table(allstraights) %>%
  as.data.frame() %>%
  filter(Freq > 2) %>%
  mutate(PID = as.character(allstraights)) %>%
  pull(PID)




#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### Complete Data #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#

completes <- df %>%
  filter(!PID %in% testPID) %>%
  filter(!PID %in% dups) %>%
  filter(!PID %in% allstraightsPID) %>%
  rename(time=`Duration (in seconds)`) %>%
  mutate(time = as.numeric(time),
         time = time/60,
         age_cat = case_when(age >= 18 & age <=24 ~ "18-24",
                             age >= 25 & age <=34 ~ "25-34",
                             age >= 35 & age <=44 ~ "35-44",
                             age >= 45 & age <=54 ~ "45-54",
                             age >= 55 & age <=64 ~ "55-64",
                             age >= 65 ~ "65+"),
         region = case_when(province == "Ontario" ~ "ON",
                            province == "British Columbia" ~ "BC",
                            province == "Quebec" ~ "QC",
                            province == "Alberta" ~ "West",
                            province == "Saskatchewan" ~ "West",
                            province == "Manitoba" ~ "West",
                            province == "New Brunswick" ~ "ATL",
                            province == "Nova Scotia" ~ "ATL",
                            province == "Newfoundland and Labrador" ~ "ATL",
                            province == "Prince Edward Island" ~ "ATL")) %>%
  mutate(across(where(is.character), ~na_if(., ""))) %>%
  mutate(ac_correct = case_when(ac == "Yellow" ~ 1,
                                TRUE ~ 0)) %>%
  drop_na(race)
  

#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### Create Scales #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
  

#### ~ Racial Resentment ~ ####

scale_rr <- completes %>%
  select(PID, contains("racialresent")) %>%
  select(-c(racialresent_DO)) %>%
  mutate(rr1 = case_when(racialresent_1 == "Strongly agree" ~ 4,
                         racialresent_1 == "Somewhat agree" ~ 3,
                         racialresent_1 == "Neither agree nor disagree" ~ 2,
                         racialresent_1 == "Somewhat disagree" ~ 1,
                         racialresent_1 == "Strongly disagree" ~ 0),
         rr2 = case_when(racialresent_2 == "Strongly agree" ~ 0,
                         racialresent_2 == "Somewhat agree" ~ 1,
                         racialresent_2 == "Neither agree nor disagree" ~ 2,
                         racialresent_2 == "Somewhat disagree" ~ 3,
                         racialresent_2 == "Strongly disagree" ~ 4),
         rr3 = case_when(racialresent_3 == "Strongly agree" ~ 0,
                         racialresent_3 == "Somewhat agree" ~ 1,
                         racialresent_3 == "Neither agree nor disagree" ~ 2,
                         racialresent_3 == "Somewhat disagree" ~ 3,
                         racialresent_3 == "Strongly disagree" ~ 4),
         rr4 = case_when(racialresent_4 == "Strongly agree" ~ 4,
                         racialresent_4 == "Somewhat agree" ~ 3,
                         racialresent_4 == "Neither agree nor disagree" ~ 2,
                         racialresent_4 == "Somewhat disagree" ~ 1,
                         racialresent_4 == "Strongly disagree" ~ 0),
         rr5 = case_when(racialresent_5 == "Strongly agree" ~ 0,
                         racialresent_5 == "Somewhat agree" ~ 1,
                         racialresent_5 == "Neither agree nor disagree" ~ 2,
                         racialresent_5 == "Somewhat disagree" ~ 3,
                         racialresent_5 == "Strongly disagree" ~ 4),
         rr6 = case_when(racialresent_6 == "Strongly agree" ~ 4,
                         racialresent_6 == "Somewhat agree" ~ 3,
                         racialresent_6 == "Neither agree nor disagree" ~ 2,
                         racialresent_6 == "Somewhat disagree" ~ 1,
                         racialresent_6 == "Strongly disagree" ~ 0)) %>%
  select(PID, contains("rr"))

summary(scale_rr)

scale_rr_eval <- select(scale_rr, -PID)

psych::alpha(scale_rr_eval)
# alpha = 0.85

scale_rr_eval$rrscale <- apply(scale_rr_eval, 1, function(x) mean(x, na.rm = TRUE))
scale_rr_eval$rrscale01 <- modEvA::range01(scale_rr_eval$rrscale, na.rm = TRUE)

quantile(scale_rr_eval$rrscale01, probs = c(0.33, 0.66), na.rm = TRUE)

plot(hist(scale_rr_eval$rrscale))
plot(hist(scale_rr_eval$rrscale01))

summary(scale_rr_eval)

rrscale <- cbind(scale_rr, scale_rr_eval) %>%
  select(PID, rrscale, rrscale01) %>%
  mutate(rrscale01_binary = case_when(rrscale01 <= 0.4583 ~ "LOW",
                                      rrscale01 > 0.4583 ~ "HIGH"),
         rrscale01_3bins = case_when(rrscale01 <= 0.333 ~ "BOTTOM",
                                     rrscale01 > 0.333 & rrscale01 <= 0.5416 ~ "MIDDLE",
                                     rrscale01 > 0.5416 ~ "TOP"))
  
head(rrscale)

#### ~ Vulnerability ~ ####


scale_vuln <- completes %>%
  select(PID, contains("empfin")) %>%
  select(-empfin_DO) %>%
  mutate(vuln1 = case_when(empfin_1 == "Strongly agree" ~ 4,
                           empfin_1 == "Somewhat agree" ~ 3,
                           empfin_1 == "Neither agree nor disagree" ~ 2,
                           empfin_1 == "Somewhat disagree" ~ 1,
                           empfin_1 == "Strongly disagree" ~ 0),
         vuln2 = case_when(empfin_2 == "Strongly agree" ~ 4,
                           empfin_2 == "Somewhat agree" ~ 3,
                           empfin_2 == "Neither agree nor disagree" ~ 2,
                           empfin_2 == "Somewhat disagree" ~ 1,
                           empfin_2 == "Strongly disagree" ~ 0),
         vuln3 = case_when(empfin_3 == "Strongly agree" ~ 0,
                           empfin_3 == "Somewhat agree" ~ 1,
                           empfin_3 == "Neither agree nor disagree" ~ 2,
                           empfin_3 == "Somewhat disagree" ~ 3,
                           empfin_3 == "Strongly disagree" ~ 4),
         vuln4 = case_when(empfin_4 == "Strongly agree" ~ 0,
                           empfin_4 == "Somewhat agree" ~ 1,
                           empfin_4 == "Neither agree nor disagree" ~ 2,
                           empfin_4 == "Somewhat disagree" ~ 3,
                           empfin_4 == "Strongly disagree" ~ 4)) %>%
  select(PID, contains("vuln"))

head(scale_vuln)

summary(scale_vuln)

scale_vuln_eval <- select(scale_vuln, -PID)

psych::alpha(scale_vuln_eval)
# alpha = 0.7


scale_vuln_eval$vulnscale <- apply(scale_vuln_eval, 1, function(x) mean(x, na.rm = TRUE))
scale_vuln_eval$vulnscale01 <- range01(scale_vuln_eval$vulnscale, na.rm = TRUE)


plot(hist(scale_vuln_eval$vulnscale))
plot(hist(scale_vuln_eval$vulnscale01))

summary(scale_vuln_eval)

vulnscale <- cbind(scale_vuln, scale_vuln_eval) %>%
  select(PID, vulnscale, vulnscale01) %>%
  mutate(vulnscale01_binary = case_when(vulnscale01 <= 0.5625 ~ "LOW",
                                        vulnscale01 > 0.5625 ~ "HIGH"))

#### ~ News Consumption ~ ####

scale_news <- completes %>%
  select(PID, contains("news_")) %>%
  select(-c(news_DO)) %>%
  pivot_longer(cols = -c(PID)) %>%
  mutate(value = case_when(value == "Never" ~ 0,
                           value == "A few times a month" ~ 1,
                           value == "A few times a week" ~ 2,
                           value == "Daily" ~ 3)) %>%
  pivot_wider(names_from = name, values_from = value)

head(scale_news)
summary(scale_news)

scale_news_eval <- select(scale_news, -PID)

psych::alpha(scale_news_eval)
# alpha = 0.66

scale_news_eval$newsscale <- apply(scale_news_eval, 1, function(x) mean(x, na.rm = TRUE))
scale_news_eval$newsscale01 <- range01(scale_news_eval$newsscale, na.rm = TRUE)


plot(hist(scale_news_eval$newsscale))
plot(hist(scale_news_eval$newsscale01))

summary(scale_news_eval)

newsscale <- cbind(scale_news, scale_news_eval) %>%
  select(PID, newsscale, newsscale01) %>%
  mutate(newsscale01_binary = case_when(newsscale01 <= 0.4762 ~ "LOW",
                                        newsscale01 > 0.4762 ~ "HIGH"))



#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### ~~~ Data for Analysis ~~~ #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#


df_analyis <- completes %>%
  left_join(rrscale, by = "PID") %>%
  left_join(vulnscale, by = "PID") %>%
  left_join(newsscale, by = "PID") %>%
  mutate(
         # sex
         sex = case_when(gender == "A man" ~ "Male",
                         TRUE ~ "Female/Other"),
         # race
         race = case_when(race == "White" ~ "White",
                          TRUE ~ "Racialized"),
         #ps degree
         psdegree = case_when(edu == "Completed CEGEP, college, or university" |
                                edu == "Graduate degree or Professional degree (e.g., MBA, JD, LLB)" ~ "High",
                              edu == "Some CEGEP, college, or university" ~ "Medium",
                              edu == "High school graduate" | edu == "Less than high school" ~ "Low"),
         psdegree = fct_relevel(psdegree, "Low", "Medium", "High"),
         # income
         income3 = case_when(income == "Less than $29,999" | income == "$30,000 to $89,999" ~ "Low",
                            income == "$90,000 to $149,999" ~ "Medium",
                            income == "$150,000 to $199,999" | income == "$200,000 to $299,999" | 
                              income == "$300,000 or more" ~ "High"),
         income3 = na_if(income3, "-99"),
         income3 = fct_relevel(income3, "Low", "Medium", "High"),
         # political interest
         polinterest = as.numeric(polinterest_1),
         # ideology
         ideo = as.numeric(ideo_1),
          #pid
         party = case_when(pid_roc == "Conservative" | pid_qc == "Conservative" ~ "CPC",
                           pid_roc == "Liberal" | pid_qc == "Liberal" ~ "LPC",
                           pid_roc == "NDP" | pid_qc == "NDP" ~ "NDP",
                           pid_roc == "Green" | pid_qc == "Green" ~ "Green",
                           pid_qc == "Bloc Québécois" ~ "BQ",
                           TRUE ~ "Other"),
         # attention check
         ac = case_when(ac == "Yellow" ~ "Pass",
                        TRUE ~ "Fail"),
         # language
         userlang = case_when(UserLanguage == "EN" ~ "English",
                              UserLanguage == "FR-CA" ~ "French"),
         # occupation
         occ = case_when(occupation == "Agriculture, fisheries, and forestry" ~ "Blue Collar",
                         occupation == "Construction" ~ "Blue Collar",
                         occupation == "Manufacturing" ~ "Blue Collar",
                         occupation == "Trades and transport" ~ "Blue Collar",
                         occupation == "Wholesale and retail" ~ "Blue Collar",
                         TRUE ~ "White Collar")
         ) %>%
  # union
  mutate(union = na_if(union, "-99")) %>%
  # immigrant
  mutate(immigrant = na_if(immigrant, "-99")) %>%
  # marital status
  mutate(maritalstatus = na_if(maritalstatus, "-99")) %>%
  # egocentric
  mutate(egocentric = na_if(egocentric, "-99"),
         egocentric = fct_relevel(egocentric, 
                                  "Got a lot worse",
                                  "Got a little worse",
                                  "Stayed about the same",
                                  "Got a little better",
                                  "Got a lot better")) %>%
  # sociotropic
  mutate(sociotropic = na_if(sociotropic, "-99"),
         sociotropic = fct_relevel(sociotropic, 
                                   "Got a lot worse",
                                   "Got a little worse",
                                   "Stayed about the same",
                                   "Got a little better",
                                   "Got a lot better")) %>%
  mutate(polinterest = na_if(polinterest, -99)) %>%
  mutate(ideo = na_if(ideo, -99))


write.csv(df_analyis, file = here("Data", "Data for Analysis.csv"))


#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### ~~~ Descriptive Tables ~~~ #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#


dtab <- df_analyis %>%
  select(sex, race, age_cat, income3, immigrant, psdegree, maritalstatus, union, occ, polinterest, party, ideo,
         rrscale01, rrscale01_binary, newsscale01, newsscale01_binary, vulnscale01, vulnscale01_binary,
         egocentric, sociotropic, region, ac, userlang)

head(dtab)

table1::label(dtab$sex) <- "Sex"
table1::label(dtab$race) <- "Race"
table1::label(dtab$age_cat) <- "Age"
table1::label(dtab$income3) <- "Income"
table1::label(dtab$immigrant) <- "Parent Immigrant"
table1::label(dtab$psdegree) <- "Education"
table1::label(dtab$maritalstatus) <- "Marital Status"
table1::label(dtab$union) <- "Union Member"
table1::label(dtab$occ) <- "Occupation Type"
table1::label(dtab$polinterest) <- "Political Interest"
table1::label(dtab$party) <- "Party ID"
table1::label(dtab$ideo) <- "Ideological Self-Placement"
table1::label(dtab$rrscale01) <- "Racial Resentment Index"
table1::label(dtab$rrscale01_binary) <- "Racial Resentment Index (Binary)"
table1::label(dtab$newsscale01) <- "Media Consumption Index"
table1::label(dtab$newsscale01_binary) <- "Media Consumption Index (Binary)"
table1::label(dtab$vulnscale01) <- "Vulnerability Index"
table1::label(dtab$vulnscale01_binary) <- "Vulnerability Index (Binary)"
table1::label(dtab$egocentric) <- "Egocentric"
table1::label(dtab$sociotropic) <- "Sociotropic"
table1::label(dtab$region) <- "Region"
table1::label(dtab$ac) <- "Attention Check"
table1::label(dtab$userlang) <- "Survey Language"

dtab_covars <- table1::table1(~sex + race + age_cat + income3 + immigrant + psdegree + maritalstatus + union + occ +
                 polinterest + party + ideo + egocentric + sociotropic + region + ac + userlang, data = dtab)

dtab_index <- table1::table1(~rrscale01 + rrscale01_binary +
                 newsscale01 + newsscale01_binary +
                 vulnscale01 + vulnscale01_binary, data = dtab)

# Convert to flextable
t1flex(dtab_covars) %>% 
  save_as_docx(path=here("Tables", "SI_Descriptive Table (Covars).docx"))

# Convert to flextable
t1flex(dtab_index) %>% 
  save_as_docx(path=here("Tables", "SI_Descriptive Table (Indices).docx"))

#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### ~~~ ANALYSIS ~~~ #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
  
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### 1. Multidimensional Preference Scaling #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#

mdpref_ind <- df_analyis %>%
  filter(ac == "Pass") %>%
  select(PID, ac, sex, race, age_cat, psdegree, party, income3, union, occ, rrscale01, newsscale01, vulnscale01, region,
         manufacturing=indpref_0_1_RANK,
         healthcare=indpref_0_2_RANK,
         agriculture=indpref_0_3_RANK,
         education=indpref_0_4_RANK,
         natresources=indpref_0_5_RANK,
         retail=indpref_0_6_RANK) %>%
  mutate(
         manufacturing = as.numeric(manufacturing),
         healthcare = as.numeric(healthcare),
         agriculture = as.numeric(agriculture),
         education = as.numeric(education),
         natresources = as.numeric(natresources),
         retail = as.numeric(retail)) %>%
  drop_na(manufacturing)

head(mdpref_ind)

summary(mdpref_ind)

indscale <- mdpref_ind %>%
  select(manufacturing, healthcare, agriculture, education, natresources, retail)


summary(indscale)


#### ~~~ Summary Table ~~~ ####

tab_indscale <- mdpref_ind %>%
  select(PID, manufacturing, healthcare, agriculture, education, natresources, retail) %>%
  pivot_longer(cols = -c(PID)) %>%
  group_by(name, value) %>%
  tally() %>%
  mutate(ntot = 1899) %>%
  mutate(npct = round(n/ntot*100, 3))

tab_indscale %>%
  filter(name == "retail")






#define columns to reverse code
reverse_cols = c("manufacturing", 
                 "healthcare", 
                 "agriculture", 
                 "education", 
                 "natresources", 
                 "retail")

#reverse code Q2 and Q5 columns
indscale[ , reverse_cols] = 6 - indscale[ , reverse_cols]

set.seed(515)
nonmetric <- mdpref(indscale, meas.lev = 2, meas.proc = 1, dimens = 2)

nonmetric$fitvector
#0.6025154 0.8851742 0.9310872 0.9665105 1.0000000 1.0000000

# check plot

meanvec <- as.matrix(apply(nonmetric$row.coords, FUN = "mean", MARGIN = 2))

rownames(nonmetric$col.coords) <- c("Manufacturing", "Healthcare", "Agriculture",
                                    "Education", "Nat. Resources", "Retail")

set.seed(515)
xyplot(V2 ~ V1,
       data = nonmetric$row.coords,
       xlab = "",
       ylab = "",
       scales = list(draw = FALSE),
       aspect = 1,
       panel = function (x, y) {
         panel.xyplot(jitter(x, amount=0.05), 
                      jitter(y, amount=0.05), col = TRUE)
         panel.xyplot(nonmetric$col.coords$V1, nonmetric$col.coords$V2, 
                      cex = .5, col = "black", pch = 19)
         panel.text(nonmetric$col.coords$V1, nonmetric$col.coords$V2, 
                    labels = rownames(nonmetric$col.coords), pos = 4, cex = 1)
         panel.arrows(0, 0, meanvec[1], meanvec[2], col = "purple")
         panel.text(-0.1, 0.15, "Avg. Respondent", pos = 4, cex = 1, col = "purple")
         
       }
)


#### ~~ Save Circle Figure ~~ ####

png(filename = here("Figures", "MDPREF Output.png"), width = 480, height = 480, units = "px", pointsize = 12)

set.seed(515)
xyplot(V2 ~ V1,
       data = nonmetric$row.coords,
       xlab = "",
       ylab = "",
       scales = list(draw = FALSE),
       aspect = 1,
       panel = function (x, y) {
         panel.xyplot(jitter(x, amount=0.05), 
                      jitter(y, amount=0.05), col = TRUE)
         panel.xyplot(nonmetric$col.coords$V1, nonmetric$col.coords$V2, 
                      cex = .5, col = "black", pch = 19)
         panel.text(nonmetric$col.coords$V1, nonmetric$col.coords$V2, 
                    labels = rownames(nonmetric$col.coords), pos = 4, cex = 1)
         panel.arrows(0, 0, meanvec[1], meanvec[2], col = "purple")
         panel.text(-0.1, 0.15, "Avg. Respondent", pos = 4, cex = 1, col = "purple")
         
       }
)

dev.off()

#### ~ Combine Output ~ ####

mdpref_combdata <- cbind(mdpref_ind$sex, 
                         mdpref_ind$race,
                         mdpref_ind$age_cat,
                         mdpref_ind$psdegree,
                         mdpref_ind$party, 
                         mdpref_ind$income,
                         mdpref_ind$union,
                         mdpref_ind$occ,
                         mdpref_ind$rrscale01,
                         mdpref_ind$newsscale01,
                         mdpref_ind$vulnscale01,
                         mdpref_ind$region, 
                         nonmetric$row.coords)

head(mdpref_combdata)
names(mdpref_combdata) <- c("sex", "race","age", "psdegree", "party", 
                            "income", "union", "occupation", "rr", "news", "vuln", "region",
                            "V1", "V2")

head(mdpref_combdata)

mdpref_combdata <- na.omit(mdpref_combdata)

mdpref_combdata <- dummy_cols(mdpref_combdata, select_columns = "age", remove_first_dummy = TRUE)

mdpref_combdata <- mdpref_combdata %>%
  mutate(region = fct_relevel(region, "West"))

mdpref_combdata <- dummy_cols(mdpref_combdata, select_columns = "region", remove_first_dummy = TRUE)

mdpref_combdata <- dummy_cols(mdpref_combdata, select_columns = "income", remove_first_dummy = TRUE)

mdpref_combdata <- dummy_cols(mdpref_combdata, select_columns = "psdegree", remove_first_dummy = TRUE)

mdpref_combdata <- mdpref_combdata %>%
  mutate(party = fct_relevel(party, "LPC"))

mdpref_combdata <- dummy_cols(mdpref_combdata, select_columns = "party", remove_first_dummy = TRUE)

mdpref_combdata <- mdpref_combdata %>%
  mutate(sex = case_when(sex == "Male" ~ 1,
                         TRUE ~ 0),
         race = case_when(race == "White" ~ 1,
                          TRUE ~ 0),
         union = case_when(union == "Yes" ~ 1,
                           union == "No" ~ 0),
         occupation = case_when(occupation == "Blue Collar" ~ 1,
                                occupation == "White Collar" ~ 0))
  
  
  
head(mdpref_combdata)

mdpref_combdata <- mdpref_combdata %>%
  janitor::clean_names()

#### ~ Circular Regression ~ ####


###
###   Calculate angle of each row vector, in radians,
###   using 3:00 position as origin, going counter-
###   clockwise. Also declare the angle to be a 
###   circular object, and create new version of
###   data frame that omits observations with
###   missing values
###

mdpref_combdata$ang <- atan(mdpref_combdata$v2/mdpref_combdata$v1) + ((mdpref_combdata$v1 < 0) * pi) +
  ((mdpref_combdata$v2 <0 & mdpref_combdata$v1 > 0) * 2 * pi)

mdpref_combdata$angle <- circular(mdpref_combdata$ang)

circreg2 <- na.omit(mdpref_combdata)

###
###   Center continuous variables
###

circreg2$rr <- circreg2$rr - mean(circreg2$rr)
circreg2$news <- circreg2$news - mean(circreg2$news)
circreg2$vuln <- circreg2$vuln - mean(circreg2$vuln)

###
###   Carry out circular regression, producing
###   results reported in Table 3
###

creg1 <- lm.circular(y = circreg2$angle,
                     x = cbind(circreg2$sex,
                               circreg2$race,
                               circreg2$age_25_34,
                               circreg2$age_35_44,
                               circreg2$age_45_54,
                               circreg2$age_55_64,
                               circreg2$age_65,
                               circreg2$income_medium,
                               circreg2$income_high,
                               circreg2$psdegree_medium,
                               circreg2$psdegree_high,
                               circreg2$party_bq,
                               circreg2$party_cpc,
                               circreg2$party_green,
                               circreg2$party_ndp,
                               circreg2$party_other,
                               circreg2$union,
                               circreg2$occupation,
                               circreg2$rr,
                               circreg2$news,
                               circreg2$vuln,
                               circreg2$region_atl,
                               circreg2$region_bc,
                               circreg2$region_on,
                               circreg2$region_qc),
                     type = "c-l",
                     init = rep(.008, 25),
                     verbose = T
)

# 470 iterations

creg1

# mu = 2.3 (approximately 10:00)


#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### 2. Experiment (Info) #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#



exp2 <- df_analyis %>%
  filter(ac == "Pass")  %>%
  mutate(expgroup2 = case_when(Experiment2Introduction_DO == "exp2_intro_withstats" ~ "Gendered Impact Information",
                               Experiment2Introduction_DO == "exp2_intro_nostats" ~ "Control"),
         exp2_gains_man = as.numeric(exp2_gains_1),
         exp2_losses_man = as.numeric(exp2_losses_1),
         exp2_gains_healthcare = as.numeric(exp2_gains_2),
         exp2_losses_healthcare = as.numeric(exp2_losses_2))

head(exp2)

summary(exp2)

#### ~ Results ~ ####

m1 <- lm(exp2_gains_healthcare ~ expgroup2, data = exp2)
summary(m1)

m1a <- lm(exp2_gains_healthcare ~ expgroup2*sex, data = exp2)
summary(m1a)
m1a_man <- lm(exp2_gains_man ~ expgroup2*sex, data = exp2)
summary(m1a_man)

m2 <- lm(exp2_losses_healthcare ~ expgroup2, data = exp2)
summary(m2)

m2a <- lm(exp2_losses_healthcare ~ expgroup2*sex, data = exp2)
summary(m2a)
m2a_man <- lm(exp2_losses_man ~ expgroup2*sex, data = exp2)
summary(m2a_man)

plot_predictions(m1a, by = c("expgroup2", "sex"), conf_level = 0.95, draw = TRUE)
plot_predictions(m2a, by = c("expgroup2", "sex"), conf_level = 0.95, draw = TRUE)

nobs(m1a)
nobs(m2a)


## Stat. sig. for plot

## job gains

avg_comparisons(
  m1a, 
  newdata = datagrid(sex="Male"),
  variables = list(expgroup2=c("Control", "Gendered Impact Information"))
)

avg_comparisons(
  m1a, 
  newdata = datagrid(sex="Female/Other"),
  variables = list(expgroup2=c("Control", "Gendered Impact Information"))
)

avg_comparisons(m1a, variables = "sex")

## job losses

avg_comparisons(
  m2a, 
  newdata = datagrid(sex="Male"),
  variables = list(expgroup2=c("Control", "Gendered Impact Information"))
)

avg_comparisons(
  m2a, 
  newdata = datagrid(sex="Female/Other"),
  variables = list(expgroup2=c("Control", "Gendered Impact Information"))
)

avg_comparisons(m2a, variables = "sex")

#### ~ Plots ~ ####

exp2_plot1 <- plot_predictions(m1a, by = c("expgroup2", "sex"), conf_level = 0.95, draw = FALSE) %>%
  mutate(sex = case_when(sex == "Male" ~ "Man",
                         sex == "Female/Other" ~ "Woman/Other"),
         sex = fct_relevel(sex, "Woman/Other"))

exp2_plot1 <- ggplot(exp2_plot1, aes(x = sex, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 5, color = "indianred", linetype="dashed", linewidth = 1.5) +
  geom_pointrange(aes(color = expgroup2, shape = expgroup2), stat = "identity", position = position_dodge2(width = .5), size = 1) +
  scale_y_continuous(limits = c(3.5, 6.5)) +      
  scale_color_manual(values = c("#d55e00", "#0072b2")) +
  geom_signif(comparisons=list(c("Woman/Other", "Man")),
              inherit.aes = TRUE,
              map_signif_level = FALSE,
              annotations="***",
              y_position = 6.2, 
              tip_length = 0.1, 
              vjust=0,
              size = 1.5) +
  geom_signif(
    y_position = c(5.85, 5.4), xmin = c(0.8, 1.8), xmax = c(1.2, 2.2),
    annotation = c("NS", "NS"), tip_length = 0
  ) +
  labs(
    x = "Respondent Gender",
    y = "Desired Healthcare Job Gains (0-10)",
    title=bquote("Healthcare Job"~bold("Gains"))
  ) +
  theme_minimal() +
  theme(
    legend.title = element_blank()
  )

exp2_plot1


exp2_plot2 <- plot_predictions(m2a, by = c("expgroup2", "sex"), conf_level = 0.95, draw = FALSE) %>%
  mutate(sex = case_when(sex == "Male" ~ "Man",
                         sex == "Female/Other" ~ "Woman/Other"),
         sex = fct_relevel(sex, "Woman/Other"))

exp2_plot2 <- ggplot(exp2_plot2, aes(x = sex, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 5, color = "indianred", linetype="dashed", linewidth = 1.5) +
  geom_pointrange(aes(color = expgroup2, shape = expgroup2), stat = "identity", position = position_dodge2(width = .5), size = 1) +
  scale_y_continuous(limits = c(3.5, 6.5)) +      
  scale_color_manual(values = c("#d55e00", "#0072b2")) +
  geom_signif(comparisons=list(c("Woman/Other", "Man")),
              inherit.aes = TRUE,
              map_signif_level = FALSE,
              annotations="***",
              y_position = 4.75, 
              tip_length = 0.1, 
              vjust=0,
              size = 1.5) +
  geom_signif(
    y_position = c(4.25, 4.5), xmin = c(0.8, 1.8), xmax = c(1.2, 2.2),
    annotation = c("NS", "NS"), tip_length = 0
  ) +
  labs(
    x = "Respondent Gender",
    y = "Desired Healthcare Job Losses (0-10)",
    #title = "Desired Healthcare Job Losses (0-10)"
    title=bquote("Healthcare Job"~bold("Losses"))
  ) +
  theme_minimal() +
  theme(
    legend.title = element_blank()
  )

exp2_plot2

combo1 <- ggpubr::ggarrange(exp2_plot1, exp2_plot2,
                            nrow = 1,
                            ncol = 2,
                            common.legend = TRUE,
                            legend = "bottom")


combo1

ggsave(combo1, filename = here("Figures", "Gender, Race, and Trade Survey - Experiment 2.png"), 
       scale = 1.2, dpi = 600, bg="white")


#### ~ Conditional Effects on Occupation Type ~ ####
m3a <- lm(exp2_gains_healthcare ~ expgroup2*sex*occ, data = exp2)
m4a <- lm(exp2_losses_healthcare ~ expgroup2*sex*occ, data = exp2)

summary(m3a)
summary(m4a)


plot_predictions(m3a, by = c("sex", "expgroup2", "occ"), conf_level = 0.95, draw = TRUE)
plot_predictions(m4a, by = c("sex", "expgroup2", "occ"), conf_level = 0.95, draw = TRUE)





## Stat. sig. for plot

## job gains

avg_comparisons(
  m3a, 
  newdata = datagrid(occ="White Collar"),
  variables = list(expgroup2=c("Control", "Gendered Impact Information"),
                   sex=c("Male","Female/Other"))
)

avg_comparisons(
  m3a, 
  newdata = datagrid(occ="Blue Collar"),
  variables = list(expgroup2=c("Control", "Gendered Impact Information"),
                   sex=c("Male","Female/Other"))
)

## job losses

avg_comparisons(
  m4a, 
  newdata = datagrid(occ="White Collar"),
  variables = list(expgroup2=c("Control", "Gendered Impact Information"),
                   sex=c("Male","Female/Other"))
)

avg_comparisons(
  m4a, 
  newdata = datagrid(occ="Blue Collar"),
  variables = list(expgroup2=c("Control", "Gendered Impact Information"),
                   sex=c("Male","Female/Other"))
)


#### ~ Plots ~ ####

# gains (BC)

exp2_plot1alt <- plot_predictions(m3a, by = c("expgroup2", "sex", "occ"), conf_level = 0.95, draw = FALSE) %>%
  mutate(sex = case_when(sex == "Male" ~ "Man",
                         sex == "Female/Other" ~ "Woman/Other"),
         sex = fct_relevel(sex, "Woman/Other"))

exp2_plot1alt_bc <- ggplot(filter(exp2_plot1alt, occ == "Blue Collar"), aes(x = sex, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 5, color = "indianred", linetype="dashed", linewidth = 1.5) +
  geom_pointrange(aes(color = expgroup2, shape = expgroup2), stat = "identity", position = position_dodge2(width = .5), size = 1) +
  scale_y_continuous(limits = c(3.5, 6.5)) +      
  scale_color_manual(values = c("#d55e00", "#0072b2")) +
  geom_signif(comparisons=list(c("Woman/Other", "Man")),
              inherit.aes = TRUE,
              map_signif_level = FALSE,
              annotations=c("NS (p = 0.102)"),
              y_position = 6.4, 
              tip_length = 0.1, 
              vjust=0,
              size = 1.5) +
  geom_signif(
    y_position = c(6.15, 5.6), xmin = c(0.8, 1.8), xmax = c(1.2, 2.2),
    annotation = c("NS", "NS"), tip_length = 0
  ) +
  labs(
    x = "Respondent Gender",
    y = "Desired Healthcare Job Gains (0-10)",
    title=bquote("Healthcare Job"~bold("Gains")~"(Blue Collar)")
  ) +
  theme_bw() +
  theme(
    legend.title = element_blank()
  )

exp2_plot1alt_bc

# gains (WC)

exp2_plot1alt_wc <- ggplot(filter(exp2_plot1alt, occ == "White Collar"), aes(x = sex, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 5, color = "indianred", linetype="dashed", linewidth = 1.5) +
  geom_pointrange(aes(color = expgroup2, shape = expgroup2), stat = "identity", position = position_dodge2(width = .5), size = 1) +
  scale_y_continuous(limits = c(3.5, 6.5)) +      
  scale_color_manual(values = c("#d55e00", "#0072b2")) +
  geom_signif(comparisons=list(c("Woman/Other", "Man")),
              inherit.aes = TRUE,
              map_signif_level = FALSE,
              annotations=c("***"),
              y_position = 6.4, 
              tip_length = 0.1, 
              vjust=0,
              size = 1.5) +
  geom_signif(
    y_position = c(5.9, 5.4), xmin = c(0.8, 1.8), xmax = c(1.2, 2.2),
    annotation = c("NS", "NS"), tip_length = 0
  ) +
  labs(
    x = "Respondent Gender",
    y = "Desired Healthcare Job Gains (0-10)",
    title=bquote("Healthcare Job"~bold("Gains")~"(White Collar)")
  ) +
  theme_bw() +
  theme(
    legend.title = element_blank()
  )

exp2_plot1alt_wc

###

exp2_plot2alt <- plot_predictions(m4a, by = c("expgroup2", "sex", "occ"), conf_level = 0.95, draw = FALSE) %>%
  mutate(sex = case_when(sex == "Male" ~ "Man",
                         sex == "Female/Other" ~ "Woman/Other"),
         sex = fct_relevel(sex, "Woman/Other"))

# losses (BC)

exp2_plot2alt_bc <- ggplot(filter(exp2_plot2alt, occ == "Blue Collar"), aes(x = sex, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 5, color = "indianred", linetype="dashed", linewidth = 1.5) +
  geom_pointrange(aes(color = expgroup2, shape = expgroup2), stat = "identity", position = position_dodge2(width = .5), size = 1) +
  scale_y_continuous(limits = c(3, 6.5)) +      
  scale_color_manual(values = c("#d55e00", "#0072b2")) +
  geom_signif(comparisons=list(c("Woman/Other", "Man")),
              inherit.aes = TRUE,
              map_signif_level = FALSE,
              annotations=c("NS"),
              y_position = 5.15, 
              tip_length = 0.1, 
              vjust=0,
              size = 1.5) +
  geom_signif(
    y_position = c(4.5, 4.8), xmin = c(0.8, 1.8), xmax = c(1.2, 2.2),
    annotation = c("NS", "NS"), tip_length = 0
  ) +
  labs(
    x = "Respondent Gender",
    y = "Desired Healthcare Job Gains (0-10)",
    title=bquote("Healthcare Job"~bold("Losses")~"(Blue Collar)")
  ) +
  theme_bw() +
  theme(
    legend.title = element_blank()
  )

exp2_plot2alt_bc

# losses (WC)

exp2_plot2alt_wc <- ggplot(filter(exp2_plot2alt, occ == "White Collar"), aes(x = sex, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_hline(yintercept = 5, color = "indianred", linetype="dashed", linewidth = 1.5) +
  geom_pointrange(aes(color = expgroup2, shape = expgroup2), stat = "identity", position = position_dodge2(width = .5), size = 1) +
  scale_y_continuous(limits = c(3.5, 6.5)) +      
  scale_color_manual(values = c("#d55e00", "#0072b2")) +
  geom_signif(comparisons=list(c("Woman/Other", "Man")),
              inherit.aes = TRUE,
              map_signif_level = FALSE,
              annotations=c("*"),
              y_position = 5.15, 
              tip_length = 0.1, 
              vjust=0,
              size = 1.5) +
  geom_signif(
    y_position = c(4.3, 4.6), xmin = c(0.8, 1.8), xmax = c(1.2, 2.2),
    annotation = c("NS", "NS"), tip_length = 0
  ) +
  labs(
    x = "Respondent Gender",
    y = "Desired Healthcare Job Gains (0-10)",
    title=bquote("Healthcare Job"~bold("Losses")~"(White Collar)")
  ) +
  theme_bw() +
  theme(
    legend.title = element_blank()
  )

exp2_plot2alt_wc


### combine

combo2 <- ggpubr::ggarrange(exp2_plot1alt_wc, exp2_plot2alt_wc, exp2_plot1alt_bc, exp2_plot2alt_bc,
                            nrow = 2,
                            ncol = 2,
                            common.legend = TRUE,
                            legend = "bottom")


combo2

ggsave(combo2, filename = here("Figures", "Gender, Race, and Trade Survey - Experiment 2 (Conditional on Occupation Type).png"),
       scale = 1.8, dpi = 600, bg="white")



#### ~ Reg Table ~ ####

modlist_exp2 <- list(
  'Healthcare Gains' = m1a, 
  'Manufacuturing Gains' = m1a_man, 
  'Healthcare Losses' = m2a,
  'Manufacuturing Losses' = m2a_man)

modelsummary(models = modlist_exp2,
             stars = TRUE)

cm <- c('expgroup2Gendered Impact Information'    = 'Treatment (Gendered Impact Info.)',
        'sexMale'    = 'Sex (Male)',
        'expgroup2Gendered Impact Information:sexMale' = 'Treatment x Sex',
        '(Intercept)' = 'Constant')

modelsummary(models = modlist_exp2,
             output = "default",
             stars = TRUE,
             estimate = "{estimate} ({std.error}){stars}",
             statistic = NULL,
             coef_map = cm,
             gof_map = c("nobs", "r.squared"))

# save output

modelsummary(models = modlist_exp2,
             output = here("Tables", "Gendered Impact Information Experiment Regression Results.docx"),
             stars = TRUE,
             estimate = "{estimate} ({std.error}){stars}",
             statistic = NULL,
             coef_map = cm,
             gof_map = c("nobs", "r.squared"))


#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### Experiment Prep ####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#


expgroups <- df_analyis %>%
  filter(ac == "Pass")  %>%
  mutate(exp1group = case_when(FL_28_DO == "FL_29" ~ "Group 1",
                               FL_28_DO == "FL_30" ~ "Group 2",
                               FL_28_DO == "FL_31" ~ "Group 3",
                               FL_28_DO == "FL_32" ~ "Group 4")) %>%
  select(PID, exp1group, contains("exp_group")) %>%
  pivot_longer(cols = -c(PID, exp1group)) %>%
  mutate(value = na_if(value, "")) %>%
  mutate(exp1check = case_when(exp1group == "Group 1" & value =="Construction" |
                                 exp1group == "Group 1" & value =="Natural Resource Extraction" |
                                 exp1group == "Group 1" & value =="Manufacturing" |
                                 exp1group == "Group 2" & value =="Construction" |
                                 exp1group == "Group 2" & value =="Natural Resource Extraction" |
                                 exp1group == "Group 2" & value =="Manufacturing" |
                                 exp1group == "Group 3" & value =="Hospitality" |
                                 exp1group == "Group 3" & value =="Retail Sales" |
                                 exp1group == "Group 3" & value =="Healthcare" |
                                 exp1group == "Group 4" & value =="Hospitality" |
                                 exp1group == "Group 4" & value =="Retail Sales" |
                                 exp1group == "Group 4" & value =="Healthcare" ~ 1,
                               TRUE ~ 0)) %>%
  drop_na(value) %>%
  group_by(PID, exp1group) %>%
  summarise(exp1check = mean(exp1check)) %>%
  mutate(exp1treat_race = case_when(exp1group == "Group 2" | exp1group == "Group 4" ~ "Racialized",
                                    exp1group == "Group 1" | exp1group == "Group 3" ~ "White"),
         exp1treat_gender = case_when(exp1group == "Group 1" | exp1group == "Group 2" ~ "Masculine",
                                      exp1group == "Group 3" | exp1group == "Group 4" ~ "Feminine"))

head(expgroups)



#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### 3. Experiment (Trade Views) #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#


tradeviews <- df_analyis %>%
  filter(ac == "Pass")  %>%
  left_join(expgroups, by = "PID") %>%
  mutate(tradeviews_promotetrade = case_when(exp1_mutz1 == "Strongly encouraged" ~ 2,
                                             exp1_mutz1 == "Somewhat encouraged" ~ 1,
                                             exp1_mutz1 == "Neither encouraged nor discouraged" ~ 0,
                                             exp1_mutz1 == "Somewhat discouraged" ~ -1,
                                             exp1_mutz1 == "Strongly discouraged" ~ -2),
         #tradeviews_promotetrade = na_if(tradeviews_promotetrade, ""),
         tradeviews_promotetrade_ord = fct_relevel(exp1_mutz1, "Strongly discouraged", "Somewhat discouraged", "Neither encouraged nor discouraged",
                                             "Somewhat encouraged", "Strongly encouraged"),
         party = fct_relevel(party, "LPC")) %>%
  drop_na(tradeviews_promotetrade)



#### ~ Results ~ ####

mod_tradeviews <- lm(tradeviews_promotetrade ~ exp1treat_race*exp1treat_gender + 
                 sex + age_cat + income3 + race + psdegree + union + party + vulnscale01 + rrscale01 + newsscale01 + region, 
               data = tradeviews)

summary(mod_tradeviews)

mod_tradeviews_exp <- lm(tradeviews_promotetrade ~ exp1treat_race*exp1treat_gender, data = tradeviews)

summary(mod_tradeviews_exp)

mod_tradeviews_s <- lm(tradeviews_promotetrade ~ exp1treat_race*exp1treat_gender*sex + 
                 age_cat + income3 + race + psdegree + union + party + vulnscale01 + rrscale01 + newsscale01 + region, 
               data = tradeviews)

summary(mod_tradeviews_s)


mod_tradeviews_r <- lm(tradeviews_promotetrade ~ exp1treat_race*exp1treat_gender*race + 
                 sex + age_cat + income3 + psdegree + union + party + vulnscale01 + rrscale01 + newsscale01 + region, 
               data = tradeviews)

summary(mod_tradeviews_r)


#### ~ Plots ~ ####

df_coef <- broom::tidy(mod_tradeviews, conf.int = TRUE, conf.level = 0.95) %>%
  filter(term != "(Intercept)") %>%
  filter(!str_detect(term, "region")) %>%
  mutate(statsig = case_when(conf.low > 0 | conf.high < 0 ~ "Yes",
                             TRUE ~ "No")) %>%
  mutate(plotord = c(1, 
                     2, 
                     4,
                     6,
                     7,
                     8,
                     9,
                     10,
                     11,
                     12,
                     5,
                     13,
                     14,
                     15,
                     16,
                     17,
                     18,
                     19,
                     20,
                     21,
                     22,
                     23,
                     3)) %>%
  mutate(termname = c(
    'Exp. Race (White)',
    'Exp. Gender (Masculine)',
    'Gender (Man)',
    'Age (25-34)',
    'Age (35-44)',
    'Age (45-54)',
    'Age (55-64)',
    'Age (65+)',
    'Income (Medium)',
    'Income (High)',
    'Race (White)',
    'Education (Medium)',
    'Education (High)',
    'Union Member (Yes)',
    'Party ID (BQ)',
    'Party ID (CPC)',
    'Party ID (Green)',
    'Party ID (NDP)',
    'Party ID (Other)',
    'Financial Vulnerability (0-1)',
    'Racial Resentment (0-1)',
    'Media Consumption (0-1)',
    'Exp. Race x Exp. Gender'
  ))

coefplot <- ggplot(df_coef, aes(y = reorder(termname, plotord), x = estimate, xmin = conf.low, xmax = conf.high)) +
  geom_point(aes(alpha = statsig)) +
  geom_linerange(aes(alpha = statsig)) +
  scale_alpha_manual(values = c(0.2, 1)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "indianred") +
  theme_minimal() +
  theme(legend.position = "none") +
  labs(
    x = "Estimate",
    y = ""
  )

coefplot

ggsave(coefplot, filename = here("Figures", "Tradeviews Question Coefficient Plot.png"), 
       scale = 1.2, dpi = 600, bg="white")

#### ~ Reg Table ~ ####

modlist_tradeviews <- list(
  mod_tradeviews_exp,
  mod_tradeviews,
  mod_tradeviews_s,
  mod_tradeviews_r
)

modelsummary(models = modlist_tradeviews,
             stars = TRUE)

cm <- c('exp1treat_raceWhite' = 'Exp. Race (White)',
        'exp1treat_genderMasculine' = 'Exp. Gender (Masculine)',
        'exp1treat_raceWhite:exp1treat_genderMasculine' = 'Exp. Race x Exp. Gender',
        
        'exp1treat_raceWhite:sexMale' = 'Exp. Race x Sex (Male)',
        'exp1treat_genderMasculine:sexMale' = 'Exp. Gender x Sex (Male)',
        
        'exp1treat_raceWhite:raceWhite' = 'Exp. Race x Race (White)',
        'exp1treat_genderMasculine:raceWhite' = 'Exp. Gender x Race (White)',
        
        'exp1treat_raceWhite:exp1treat_genderMasculine:sexMale' = "Exp. Race x Exp. Gender x Sex (Male)",
        'exp1treat_raceWhite:exp1treat_genderMasculine:raceWhite' = "Exp. Race x Exp. Gender x Race (White)",
        
        'sexMale' = 'Sex (Male)',
        'raceWhite' = 'Race (White)',
        'age_cat25-34' = 'Age (25-34)',
        'age_cat35-44' = 'Age (35-44)',
        'age_cat45-54' = 'Age (45-54)',
        'age_cat55-64' = 'Age (55-64)',
        'age_cat65+' = 'Age (65+)',
        'income3Medium' = 'Income (Medium)',
        'income3High' = 'Income (High)',
        'psdegreeMedium' = 'Education (Medium)',
        'psdegreeHigh' = 'Education (High)',
        'unionYes' = 'Union Member (Yes)',
        'partyCPC' = 'Party ID (CPC)',
        'partyNDP' = 'Party ID (NDP)',
        'partyBQ' = 'Party ID (BQ)',
        'partyGreen' = 'Party ID (Green)',
        'partyOther' = 'Party ID (Other)',
        'vulnscale01' = 'Financial Vulnerability (0-1)',
        'rrscale01' = 'Racial Resentment (0-1)',
        'newsscale01' = 'Media Consumption (0-1)',
        '(Intercept)' = 'Constant')

cm_exp <- c('exp1treat_raceWhite' = 'Exp. Race (White)',
        'exp1treat_genderMasculine' = 'Exp. Gender (Masculine)',
        'exp1treat_raceWhite:exp1treat_genderMasculine' = 'Exp. Race x Exp. Gender',
        
        'exp1treat_raceWhite:sexMale' = 'Exp. Race x Sex (Male)',
        'exp1treat_genderMasculine:sexMale' = 'Exp. Gender x Sex (Male)',
        
        'exp1treat_raceWhite:raceWhite' = 'Exp. Race x Race (White)',
        'exp1treat_genderMasculine:raceWhite' = 'Exp. Gender x Race (White)',
        
        'exp1treat_raceWhite:exp1treat_genderMasculine:sexMale' = "Exp. Race x Exp. Gender x Sex (Male)",
        'exp1treat_raceWhite:exp1treat_genderMasculine:raceWhite' = "Exp. Race x Exp. Gender x Race (White)",
        
        'sexMale' = 'Sex (Male)',
        'raceWhite' = 'Race (White)')

modelsummary(models = modlist_tradeviews,
             output = "default",
             stars = TRUE,
             estimate = "{estimate} ({std.error}){stars}",
             statistic = NULL,
             coef_map = cm,
             gof_map = c("nobs", "r.squared"))

# save output

modelsummary(models = modlist_tradeviews,
             output = here("Tables", "Tradeviews Question Experiment Regression Results.docx"),
             stars = TRUE,
             estimate = "{estimate} ({std.error}){stars}",
             statistic = NULL,
             coef_map = cm,
             gof_map = c("nobs", "r.squared"))


modelsummary(models = modlist_tradeviews,
             output = "default",
             stars = TRUE,
             estimate = "{estimate} ({std.error}){stars}",
             statistic = NULL,
             coef_map = cm_exp,
             gof_map = c("nobs", "r.squared"))

# save output

modelsummary(models = modlist_tradeviews,
             output = here("Tables", "Tradeviews Question Experiment Regression Results (Exclude Covars).docx"),
             stars = TRUE,
             estimate = "{estimate} ({std.error}){stars}",
             statistic = NULL,
             coef_map = cm_exp,
             gof_map = c("nobs", "r.squared"))



#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### 4. Experiment (Classification Primes) #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#

exp1questions <- df_analyis %>%
  filter(ac == "Pass")  %>%
  select(PID, contains("exp1_mat")) %>%
  select(-c(exp1_mat2_DO, exp1_mat1_DO)) %>%
  pivot_longer(cols = -c(PID)) %>%
  mutate(value = na_if(value, "")) %>%
  mutate(value = case_when(value == "Strongly agree" ~ 2,
                           value == "Somewhat agree" ~ 1,
                           value == "Neither agree nor disagree" ~ 0,
                           value == "Somewhat disagree" ~ -1,
                           value == "Strongly disagree" ~ -2)) %>%
  mutate(name = case_when(name == "exp1_mat1_1" ~ "exp1_trade",
                          name == "exp1_mat1_2" ~ "exp1_retrain",
                          name == "exp1_mat1_3" ~ "exp1_unemp",
                          name == "exp1_mat2_1" ~ "exp1_bene_canada",
                          name == "exp1_mat2_2" ~ "exp1_bene_likeme",
                          name == "exp1_mat2_3" ~ "exp1_bene_me")) %>%
  pivot_wider(id_cols = PID, names_from = name, values_from = value)

head(exp1questions)

exp1 <- df_analyis %>%
  filter(ac == "Pass")  %>%
  left_join(expgroups, by = "PID") %>%
  left_join(exp1questions, by = "PID") %>%
  ungroup()


#### ~ Treatments Only ~ ####

mod1 <- lm(exp1_trade ~ exp1treat_race*exp1treat_gender, data = exp1)

summary(mod1)

mod2 <- lm(exp1_retrain ~ exp1treat_race*exp1treat_gender, data = exp1)

summary(mod2)

mod3 <- lm(exp1_unemp ~ exp1treat_race*exp1treat_gender, data = exp1)

summary(mod3)

mod4 <- lm(exp1_bene_canada ~ exp1treat_race*exp1treat_gender, data = exp1)

summary(mod4)

mod5 <- lm(exp1_bene_likeme ~ exp1treat_race*exp1treat_gender, data = exp1)

summary(mod5)

mod6 <- lm(exp1_bene_me ~ exp1treat_race*exp1treat_gender, data = exp1)

summary(mod6)

#### ~ X Sex ~ ####

mod1s <- lm(exp1_trade ~ exp1treat_race*exp1treat_gender*sex, data = exp1)

summary(mod1s)

mod2s <- lm(exp1_retrain ~ exp1treat_race*exp1treat_gender*sex, data = exp1)

summary(mod2s)

mod3s <- lm(exp1_unemp ~ exp1treat_race*exp1treat_gender*sex, data = exp1)

summary(mod3s)

mod4s <- lm(exp1_bene_canada ~ exp1treat_race*exp1treat_gender*sex, data = exp1)

summary(mod4s)

mod5s <- lm(exp1_bene_likeme ~ exp1treat_race*exp1treat_gender*sex, data = exp1)

summary(mod5s)

mod6s <- lm(exp1_bene_me ~ exp1treat_race*exp1treat_gender*sex, data = exp1)

summary(mod6s)

#### ~ X Race ~ ####

mod1r <- lm(exp1_trade ~ exp1treat_race*exp1treat_gender*race, data = exp1)

summary(mod1r)

mod2r <- lm(exp1_retrain ~ exp1treat_race*exp1treat_gender*race, data = exp1)

summary(mod2r)

mod3r <- lm(exp1_unemp ~ exp1treat_race*exp1treat_gender*race, data = exp1)

summary(mod3r)

mod4r <- lm(exp1_bene_canada ~ exp1treat_race*exp1treat_gender*race, data = exp1)

summary(mod4r)

mod5r <- lm(exp1_bene_likeme ~ exp1treat_race*exp1treat_gender*race, data = exp1)

summary(mod5r)

mod6r <- lm(exp1_bene_me ~ exp1treat_race*exp1treat_gender*race, data = exp1)

summary(mod6r)


#### ~ X Racial Resentment ~ ####

mod1rr <- lm(exp1_trade ~ exp1treat_race*exp1treat_gender*rrscale01, data = exp1)

summary(mod1rr)

mod2rr <- lm(exp1_retrain ~ exp1treat_race*exp1treat_gender*rrscale01, data = exp1)

summary(mod2rr)

mod3rr <- lm(exp1_unemp ~ exp1treat_race*exp1treat_gender*rrscale01, data = exp1)

summary(mod3rr)

mod4rr <- lm(exp1_bene_canada ~ exp1treat_race*exp1treat_gender*rrscale01, data = exp1)

summary(mod4rr)

mod5rr <- lm(exp1_bene_likeme ~ exp1treat_race*exp1treat_gender*rrscale01, data = exp1)

summary(mod5rr)

mod6rr <- lm(exp1_bene_me ~ exp1treat_race*exp1treat_gender*rrscale01, data = exp1)

summary(mod6rr)


#### ~ X Occupation Type ~ ####

mod1occ <- lm(exp1_trade ~ exp1treat_race*exp1treat_gender*occ, data = exp1)

summary(mod1occ)

mod2occ <- lm(exp1_retrain ~ exp1treat_race*exp1treat_gender*occ, data = exp1)

summary(mod2occ)

mod3occ <- lm(exp1_unemp ~ exp1treat_race*exp1treat_gender*occ, data = exp1)

summary(mod3occ)

mod4occ <- lm(exp1_bene_canada ~ exp1treat_race*exp1treat_gender*occ, data = exp1)

summary(mod4occ)

mod5occ <- lm(exp1_bene_likeme ~ exp1treat_race*exp1treat_gender*occ, data = exp1)

summary(mod5occ)

mod6occ <- lm(exp1_bene_me ~ exp1treat_race*exp1treat_gender*occ, data = exp1)

summary(mod6occ)



#### ~ Reg Table ~ ####

#### ~~ Treatment Only ~~ ####

modlist_exp1_treatment <- list(
  'Promote Trade' = mod1, 
  'Increase Re-Training' = mod2, 
  'Increase Unemp. Benefits' = mod3, 
  'Benefit Canada' = mod4,
  'Benefit People Like Me' = mod5,
  'Benefit Me' = mod6)

modelsummary(models = modlist_exp1_treatment,
             stars = TRUE)

cm <- c('exp1treat_raceWhite' = 'Exp. Race (White)',
        'exp1treat_genderMasculine' = 'Exp. Gender (Masculine)',
        'exp1treat_raceWhite:exp1treat_genderMasculine' = 'Exp. Race x Exp. Gender',
        '(Intercept)' = 'Constant')

modelsummary(models = modlist_exp1_treatment,
             output = "default",
             stars = TRUE,
             #estimate = "{estimate} ({std.error}){stars}",
             #statistic = NULL,
             coef_map = cm,
             gof_map = c("nobs", "r.squared"))

# save output

modelsummary(models = modlist_exp1_treatment,
             output = here("Tables", "Classification Primes Experiment Regression Results (Treatment Only).docx"),
             stars = TRUE,
             #estimate = "{estimate} ({std.error}){stars}",
             #statistic = NULL,
             coef_map = cm,
             gof_map = c("nobs", "r.squared"))


#### ~~ Sex ~~ ####

modlist_exp1_sex <- list(
  'Promote Trade' = mod1s, 
  'Increase Re-Training' = mod2s, 
  'Increase Unemp. Benefits' = mod3s, 
  'Benefit Canada' = mod4s,
  'Benefit People Like Me' = mod5s,
  'Benefit Me' = mod6s)

modelsummary(models = modlist_exp1_sex,
             stars = TRUE)

cm <- c('exp1treat_raceWhite' = 'Exp. Race (White)',
        'exp1treat_genderMasculine' = 'Exp. Gender (Masculine)',
        'sexMale' = 'Sex (Male)',
        'exp1treat_raceWhite:exp1treat_genderMasculine' = 'Exp. Race x Exp. Gender',
        'exp1treat_raceWhite:sexMale' = 'Exp. Race x Sex',
        'exp1treat_genderMasculine:sexMale' = 'Exp. Gender x Sex',
        'exp1treat_raceWhite:exp1treat_genderMasculine:sexMale' = 'Exp. Race x Exp. Gender x Sex',
        '(Intercept)' = 'Constant')

modelsummary(models = modlist_exp1_sex,
             output = "default",
             stars = TRUE,
             #estimate = "{estimate} ({std.error}){stars}",
             #statistic = NULL,
             coef_map = cm,
             gof_map = c("nobs", "r.squared"))

# save output

modelsummary(models = modlist_exp1_sex,
             output = here("Tables", "Classification Primes Experiment Regression Results (Sex).docx"),
             stars = TRUE,
             #estimate = "{estimate} ({std.error}){stars}",
             #statistic = NULL,
             coef_map = cm,
             gof_map = c("nobs", "r.squared"))


#### ~~ Race ~~ ####

modlist_exp1_race <- list(
  'Promote Trade' = mod1r, 
  'Increase Re-Training' = mod2r, 
  'Increase Unemp. Benefits' = mod3r, 
  'Benefit Canada' = mod4r,
  'Benefit People Like Me' = mod5r,
  'Benefit Me' = mod6r)

modelsummary(models = modlist_exp1_race,
             stars = TRUE)

cm <- c('exp1treat_raceWhite' = 'Exp. Race (White)',
        'exp1treat_genderMasculine' = 'Exp. Gender (Masculine)',
        'raceWhite' = 'Race (White)',
        'exp1treat_raceWhite:exp1treat_genderMasculine' = 'Exp. Race x Exp. Gender',
        'exp1treat_raceWhite:raceWhite' = 'Exp. Race x Race',
        'exp1treat_genderMasculine:raceWhite' = 'Exp. Gender x Race',
        'exp1treat_raceWhite:exp1treat_genderMasculine:raceWhite' = 'Exp. Race x Exp. Gender x Race',
        '(Intercept)' = 'Constant')

modelsummary(models = modlist_exp1_race,
             output = "default",
             stars = TRUE,
             #estimate = "{estimate} ({std.error}){stars}",
             #statistic = NULL,
             coef_map = cm,
             gof_map = c("nobs", "r.squared"))

# save output

modelsummary(models = modlist_exp1_race,
             output = here("Tables", "Classification Primes Experiment Regression Results (Race).docx"),
             stars = TRUE,
             #estimate = "{estimate} ({std.error}){stars}",
             #statistic = NULL,
             coef_map = cm,
             gof_map = c("nobs", "r.squared"))


#### ~~ Racial Resentment ~~ ####

modlist_exp1_rr <- list(
  'Promote Trade' = mod1rr, 
  'Increase Re-Training' = mod2rr, 
  'Increase Unemp. Benefits' = mod3rr, 
  'Benefit Canada' = mod4rr,
  'Benefit People Like Me' = mod5rr,
  'Benefit Me' = mod6rr)

modelsummary(models = modlist_exp1_rr,
             stars = TRUE)

cm <- c('exp1treat_raceWhite' = 'Exp. Race (White)',
        'exp1treat_genderMasculine' = 'Exp. Gender (Masculine)',
        'rrscale01' = 'Racial Resentment',
        'exp1treat_raceWhite:exp1treat_genderMasculine' = 'Exp. Race x Exp. Gender',
        'exp1treat_raceWhite:rrscale01' = 'Exp. Race x Racial Resentment',
        'exp1treat_genderMasculine:rrscale01' = 'Exp. Gender x Racial Resentment',
        'exp1treat_raceWhite:exp1treat_genderMasculine:rrscale01' = 'Exp. Race x Exp. Gender x Racial Resentment',
        '(Intercept)' = 'Constant')

modelsummary(models = modlist_exp1_rr,
             output = "default",
             stars = TRUE,
             #estimate = "{estimate} ({std.error}){stars}",
             #statistic = NULL,
             coef_map = cm,
             gof_map = c("nobs", "r.squared"))

# save output

modelsummary(models = modlist_exp1_rr,
             output = here("Tables", "Classification Primes Experiment Regression Results (Racial Resentment).docx"),
             stars = TRUE,
             #estimate = "{estimate} ({std.error}){stars}",
             #statistic = NULL,
             coef_map = cm,
             gof_map = c("nobs", "r.squared"))


#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### ~ Plots (Treatment Effects) ~ #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#

pm1 <- broom::tidy(mod1, conf.int = TRUE) %>%
  mutate(dv = "Promote Trade") %>%
  filter(term != "(Intercept)") %>%
  mutate(termpretty = case_when(term == "exp1treat_raceWhite" ~ "Worker Prime (White)",
                                term == "exp1treat_genderMasculine" ~ "Industry Prime (Masculine)",
                                term == "exp1treat_raceWhite:exp1treat_genderMasculine" ~ "Worker Prime X Industry Prime"))

pm2 <- broom::tidy(mod2, conf.int = TRUE) %>%
  mutate(dv = "Increase Re-Training") %>%
  filter(term != "(Intercept)") %>%
  mutate(termpretty = case_when(term == "exp1treat_raceWhite" ~ "Worker Prime (White)",
                                term == "exp1treat_genderMasculine" ~ "Industry Prime (Masculine)",
                                term == "exp1treat_raceWhite:exp1treat_genderMasculine" ~ "Worker Prime X Industry Prime"))

pm5 <- broom::tidy(mod5, conf.int = TRUE) %>%
  mutate(dv = "Benefit People Like Me") %>%
  filter(term != "(Intercept)") %>%
  mutate(termpretty = case_when(term == "exp1treat_raceWhite" ~ "Worker Prime (White)",
                                term == "exp1treat_genderMasculine" ~ "Industry Prime (Masculine)",
                                term == "exp1treat_raceWhite:exp1treat_genderMasculine" ~ "Worker Prime X Industry Prime"))


plotdat <- rbind(pm1, pm2, pm5)

plotexp2 <- ggplot(plotdat, aes(x = estimate, y = termpretty, xmin = conf.low, xmax = conf.high)) +
  geom_point() +
  geom_linerange() +
  geom_vline(xintercept = 0, linetype = "dashed", color = "indianred") +
  facet_wrap(~dv, ncol = 1) +
  theme_bw() +
  theme(
    axis.title.x = element_text(hjust = 0.5),
    axis.title.y = element_text(hjust = 1),
    strip.text = element_text(face = "bold",
                              size = rel(0.75), hjust = 0.5),
    strip.background = element_rect(fill = "grey90", color = NA)
  )+
  labs(
    x = "Estimate",
    y = ""
  )

plotexp2

ggsave(plotexp2, filename = here("Figures", "Effect of Image Primes (Coefficient Plot).png"), 
       scale = 1.2, dpi = 600, bg="white")

#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### ~ Plots (RR) ~ #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#

p1 <- exp1 %>%
  lm(exp1_bene_likeme ~ exp1treat_race*exp1treat_gender*rrscale01, data = .)

exp1 %>%
  filter(exp1treat_gender != "Masculine") %>%
  lm(exp1_bene_likeme ~ exp1treat_race*rrscale01, data = .) %>%
  summary()

plot_predictions(p1, condition = c("rrscale01", "exp1treat_race", "exp1treat_gender"), conf_level = 0.84, draw = TRUE)

p1plot <- plot_predictions(p1, condition = c("rrscale01", "exp1treat_race", "exp1treat_gender"), conf_level = 0.84, draw = FALSE) %>%
  select(estimate, conf.low, conf.high, rrscale01, exp1treat_race, exp1treat_gender) %>%
  mutate(outq = "Trade Benefits People \nLike Me")

p2 <- exp1 %>%
  lm(exp1_bene_canada ~ exp1treat_race*exp1treat_gender*rrscale01, data = .)

exp1 %>%
  filter(exp1treat_gender == "Masculine") %>%
  lm(exp1_bene_canada ~ exp1treat_race*rrscale01, data = .) %>%
  summary()

plot_predictions(p2, condition = c("rrscale01", "exp1treat_race", "exp1treat_gender"), conf_level = 0.84, draw = TRUE)


p2plot <- plot_predictions(p2, condition = c("rrscale01", "exp1treat_race", "exp1treat_gender"), conf_level = 0.84, draw = FALSE) %>%
  select(estimate, conf.low, conf.high, rrscale01, exp1treat_race, exp1treat_gender) %>%
  mutate(outq = "Trade Benefits Canada")



p_comb <- rbind(p1plot, p2plot) %>%
  mutate(exp1treat_race = case_when(exp1treat_race == "Racialized" ~ "Worker Prime: Racialized",
                                    exp1treat_race == "White" ~ "Worker Prime: White"),
         exp1treat_gender = case_when(exp1treat_gender == "Feminine" ~ "Industry Prime: \nFeminine",
                                      exp1treat_gender == "Masculine" ~ "Industry Prime: \nMasculine"))
head(p_comb)

p_out <- ggplot(p_comb, aes(x = rrscale01, y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_line(aes(color = exp1treat_race)) +
  geom_ribbon(aes(fill = exp1treat_race), alpha = 0.35) +
  scale_color_manual(values = c("#d55e00", "#0072b2")) +
  scale_fill_manual(values = c("#d55e00", "#0072b2")) +
  facet_grid(outq ~ exp1treat_gender) +
  theme_bw() +
  theme(
    legend.title = element_blank(),
    legend.position = "bottom",
    strip.background =element_rect(fill="#eeeeee"),
    strip.text = element_text(colour = 'black', face = "bold")
  ) +
  labs(
    x = "Racial Resentment Scale (0-1)",
    y = "Predicted Values"
  )

p_out

ggsave(p_out, filename = here("Figures", "Effect of Image Primes Conditional On Racial Resentment.png"), 
       scale = 1.2, dpi = 600, bg="white")


#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#
#### END OF SCRIPT #####
#-----------------------------------------------------------------------------------------------------------------------------#
#-----------------------------------------------------------------------------------------------------------------------------#

